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Abstract 

We come back to the analytical solution of the standard transfer problem in a stellar 
atmosphere. It consists in solving the radiative transfer equation in a homogeneous 
and isothermal plane-parallel atmosphere, with light scattering taken as isotropic 
and monochromatic. The literature on the subject is reviewed and the existing 
solution in a finite slab is improved thanks to the introduction of non classical 
auxiliary functions. Eleven- figure tables of the solution are given for typical values 
of the input parameters currently met in stellar atmospheres. 
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1 Introduction 



We consider a static and stationary plane-parallel atmosphere with isotropic 
and monochromatic scattering in every volume unit. The source function S of 
the radiation field is solution to an integral equation of the form (Ivanov [1], 
Chap. 3) 

S(r) = S*(t) + ^-^ Ei {\t- t'\)S{t') dr'. (1) 

The first term in the right-hand side describes the contribution of the (internal 
and external) sources, and the second term describes the scattering of photons. 
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In the latter, a(r) G [0, 1] is the volumic albedo for single scattering, and E\ 
is the first exponential integral function as defined by 

E 1 (r)= / exp(-r/ M )— (r > 0). (2) 
Jo u 

The optical depth variable r covers the range [0,6], where b > denotes the 
optical thickness - possibly infinite - of the atmosphere at the considered 
frequency. Frequency dependence of any quantity is not mentioned. 

Suppose that the material is in local thermodynamic equilibrium and that the 
incoming boundary conditions vanish. Then the source term S*(t) is given by 
Kirchhoff's relation 

S*(r) = [l-a(r)]B[T(r)], (3) 
where B is the Planck function for the temperature T(r). 

The purpose of the present paper is to solve analytically the problem (l)-(3) 
for constant a and T, that is in a homogeneous and isothermal atmosphere. 
Then S* — (1 — a)B{T) is a constant and the solution reads 

S(r) =QM,r)(l -o)S(T) = S(a, b, r)B(T), (4) 

where the Q- and S-functions satisfy the following integral equations: 

Q(a, b,r) = l + U b E 1 (\r- r'\)Q(a, b, r') dr', (5) 
2 Jo 

S(a, b, r) = 1 - a + - f £i(|t - r'|)S(a, b, r') dr'. (6) 
2 Jo 

Of course S — (1 — a)Q and the above problems are equivalent. In the present 
paper, we chose to solve Eq. (5) analytically and to tabulate the function S, 
since the Q-function diverges like — a when a — * 1 and b — > +oo. 

We remind the probabilistic meaning of the Q-function: Q(a, b, r) is the mean 
number of scatterings experienced by a photon emitted on the r-layer, which 
therefore has the probability P = 1 — (1 — a)Q to leave the atmosphere since 
1 — a is the photon destruction probability per scattering [2]. The function 
S = (1 — a)Q is the source function, normalized to the Planck function, of 
the radiation field propagated in a homogeneous and isothermal atmosphere 
in local thermodynamic equilibrium. 

It is important to solve Eqs. (5) or (6) very accurately for different values of 
a G [0, 1] and b > 0, in order to have a benchmark solution available to the 
validation of codes solving Eq. (1) numerically. 

The literature on the solution of problem (5) is difficult to follow; it will be seen 
further that it is incomplete in the finite case (6 < +oo). Several analytical 
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approaches exist, which all rest on a good knowledge of some classical auxiliary 
functions, firstly the functions H, X and Y and their moments. We remind 
the definition and useful properties of these functions in Section 2. Then two 
classical approaches for solving Eq. (5) are described. The former is based on 
the analytical calculation of Sobolev's resolvent function (Section 3), and 
the latter on the calculation of Ambartsumian's 5-function (Sec. 4). We will 
see in Section. 5 and 6 that it is appropriate to express the finite slab solution 
in terms of non classical auxiliary functions, which we introduce. Finally, the 
function S — (1 — a)Q is tabulated for some typical values of the parameters 
a and b in stellar atmospheres (Section 7). 



2 A reminder on the classical H-, X- and F-functions 



The analytical solution to the transfer equation involves a number of aux- 
iliary functions, defined most of time from integral equations modelling the 
multiple scattering of photons. The functions H, X and Y, whose definition 
is reminded below, are certainly the most famous functions in plane-parallel 
geometry. Their detailed study can be found in the classical monographs on ra- 
diative transfer theory: Chandrasekhar [3], Busbridge [4], . . ., Van de Hulst [5], 
Lenoble [6]. 



2. 1 Function H 



In a semi-infinite atmosphere (b = +oo), the basic auxiliary function is the 
//-function: it depends on the parameter a and on a angular variable u which 
we suppose positive. It can be defined as the unique solution, continuous over 
[0, +oo[, of either of the following equivalent equations: 



1 a r 1 dv 

— — r = 1 u H (a,v) , 

H(a,u) 2 Jo V ' J v + u' 



(7) 



or 



d dv 

T(a,u)H(a,u) — 1 + —u / H(a,v) . 

2 Jo v — u 



(8) 



These integral equations on [0, 1] define the extension of the //-function outside 
this interval. In Eq. (8), we have introduced the dispersion function 



T(a, u) 



a 

—u In 

2 



1 + u 



1 — u 



(«^±1), 



(9) 



and the integral is a Cauchy principal value when u g]0, 1[ 
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Equations (7)- (8) can be solved analytically, a possible expression for their 
solution being given below [Eq. (69)]. Putting u — > +oo into (7) and (8), and 
observing that T(a, +oo) = 1 — a, one obtains the value at infinity of the 
^/-function 

tf(a,+oo) = ^=^. (10) 
V 1 — a 

Busbridge [7] has demonstrated that the if-function is solution to the following 
integral equation: 

R(a) k(a)u a f 1 g(a,v) dv 
{,) H[a, l/k(a)\ 1 + k(a)u 2 JoH(a,v)v + u { ' 



For < a < 1, the coefficient fc(a) is the unique root in ]0, 1[ of the character- 
istic equation T (a, l/k(a)) — 0, where T(a, w) is defined by (9). We thus have 
by definition 



2k(a) 



k(a) 



1 - k(a) 



(0<a<l), (12) 



and k(0) = 1, /c(l) = by continuity. This coefficient is of great importance 
when solving the radiative transfer equation in a slab, since l/k(a) is inter- 
preted as the thermalization depth of the atmosphere (Ivanov [1], Section 3.2). 
A complete study of the /c-function can be found in Case et al. [8], where k(a) 
is denoted by Ko(c). An analytical expression will be given in Section 6, Eqs. 
(67)-(68). 

Functions R(a) and g(a,v) appearing in the right-hand side of Eq. (11) follow 
immediately from k(a) and T(a,v) since 

™=m7^i (0<o<1) - (13) 

gia ' v)= Ti(a,v)+\«/2)av ? (0<v< 1). (14) 

Equation (11) is an integral equation on [0, 1], which allows to calculate the 
if -function over [0, +oo[ once it has been solved. Putting u — > +oo on both 
sides of Eq. (11) and taking Eq. (10) into account yields 

1 1+ r R{a ) 1( N1 +-r#4d,. (15) 



v/T^ H[a,+l/k(a)\ 2 Jo H(a,v) 

It will be seen that this equation and its counterpart in a finite slab play a 
central role in the solution of Eq. (5). 
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2.2 Functions X and Y 



In a finite slab {b < +00), the function H is replaced by two functions, denoted 
by X and Y, which depend on parameters a, b and variable u > 0. These 
functions are the unique solution, continuous over ]0,+oo[, of either of the 
following couples of integral equations: 



X(a,b,u)[l - £ x (a,b,u)] + Y(a,b,u)£ Y (a,b,u) = 1, 



(16) 



X(a,b,u)£ Y (a,b,-u) + Y(a,b,u)[l — £x(a,b, —u)] = exp(—b/u) (u ^ 1), 

(17) 



or 



T(a, u)X(a, b, u) = 1 - £ x (a, b, -u) - £y (a, 6, «) exp(-6/«) (u ^ 0, 1), (18) 



T(a,w)F(a,&,w) = [l-£ x (a,b,u)]exp(-b/u)-£ Y (a,b, -u) (u^0,l). (19) 



We have introduced the functions 



£x(a,b,u) = -u X(a,b,v)- 
2 Jo v 



f y (a,&,u) = -u / y(a,6,u)- 
z Jo V 



dv 



+ u 
dv 



+ u 



(20) 



(21) 



defined everywhere except at u — — 1. The integral is calculated in the sense 
of the Cauchy principal value when u e] — 1, 0[. 

Equations (7)-(8) are retrieved by letting b — > +00 into Eqs. (16) and (18), 
since X(a, +00, w) = H(a,u), Y(a, +00, u) = 0, +00, w) = 1/H(a,u) 

and ^y(a, +oo,m) = for -u > 0. Equations (17) and (19) reduce to = 
when b — > +00. 



We introduce the moments of order of X and Y 



ao(a,b) = X(a,b,u)du, Po(a, b) — Y(a,b,u)du, 
Jo Jo 



(22) 



and let u — > +00 into Eqs. (16)-(21). Since £,x(a, b, +00) = (a/2)cto(a,b) and 
£y(a,6, +00) = (a/2)/3o(a,b), one obtains the following relation between the 
coefficients ao and (3q\ 



1 - -a (a,6) 



: 3o{a,b) 



1 2 



1 — a, 



(23) 
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and we have further, if b < +00 



X(a, b, +00) = Y(a, b, +00) 



1-a 



1 - 2^ a ° +Po)(a,b) 



1 - T^o- A)) (a, b) 



(24) 
(25) 



The generalization of Eqs. (11) and (15) in a finite atmosphere are Eqs. (34)- 
(35) and (58)- (59) of Rutily et al. [9]. We repeat here the two last-mentioned 
equations, since they will prove useful to the remainder of the present article 



1 



1 — a 



1 - -a (a,b) 

1 + R(a) {1 - £x[a, b, l/k(a)\ + £ Y [a, b, l/k(a)\ exp[-k(a)b}} 

CL 

+ 0/ 9(a,v)[l -£ x (a,b,v) + £ Y (a,b,v) exp(-b/v)]dv, (26) 
I Jo 



1 a 



l-a2 



Po(a,b) 



R(a) {[1 - &(a, b, l/k{a))\ exp(-k(a)b) + £ Y (a, b, l/k(a))} 

+ % f g(a,v){[l-Sx(a,b,v)]exp(-b/v)+S Y (a,b,v)}dv. (27) 
I Jo 



When b — > +00, Eq. (26) yields Eq. (15) again and Eq. (27) reduces to = 0. 
We have indeed (3 (a, +00) = 0, and a (a, +00) = 2/[l + y/l — a] due to Eq. 
(23). 



3 First method for solving the problem (5) 



Sobolev [10,11] has expressed the resolvent kernel of Eq. (1) in terms of the 
associated resolvent function 0, which is solution to the same equation with 
the free-term (a/2)E 1 (r), viz. 



0(a, b, r) = ^x(r) + % f E 1 (\r - r'|)0(a, b, r') dr>. 
z z Jo 



(28) 



He inferred that the Q-function can be expressed in terms of 

ip(a, b, r) = 1+ T 0(a, b, r') dr' (29) 
Jo 
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by means of the following formulae: in a semi-infinite medium [10] 

Q( a , T ) = 7j= — ^( a > r )' ( 30 ) 
V 1 — a 

and in a finite one [11] 

Q(a, b, t) = ip(a, b, b)[ip(a, b, r) + ip(a, b,b - r) - ^(a, 6, 6)], (31) 

where 

^(a, b, b) = X(a, b, +oo) = Y(a, b, +oo) (32) 
is given by Eqs. (24)-(25). 

The resolvent function has been calculated analytically by Minin [12] in a 
semi-infinite medium 

. . . k(a)R(a) r . . a f 1 g(a,v) , , , dv 
H[a,l/k(a)\ 2 Jo H{a,v) v 

and by Rogovtsov and Samson [13] in a finite slab 

0(a, b, r) = fc(a)i2(a){[l - &(a, b, l/k{a))\ exp[-k{a)r] 

-Z Y (a,b,l/k(a))exp[-k(a)(b-T)]} 

& 

+ o / 9(a,v){[l -£ x (a,b,v)]exp(-r/v) 
2 Jo 

-Z Y (a,b,v)eM-(b-T)/v]}—. (34) 



The latter expression does yield the former again when b — > +oo. The functions 
appearing in the right-hand side have been introduced previously: see Eqs. 
(12)-(14) and (20)-(21). 

Substituting these expressions into Eq. (29), one can perform the integration 
with respect to r' analytically and simplify the resulting expression of the 
■0-function with the help of the integral Eqs. (15) and (26). The result is 



ip(a, r) = 



in a semi-infinite medium, and 



R(a) a f 1 g(a, v) 

eXp[ -^ (a)T] " 2 Jo Wjhvj eM ~ T/v) dV 

(35) 



ip(a,b,r) 



1 



1 



^a {a,b) 



-R(a){[l-£x(a, b, l/k(a))\ exp[-fc(a)r] +^r(a, 6, l/k(a)) exp[-k(a)(b-r)]} 

a rl 

- 2 J 9(a, v){[l- £x(a, b, v)] exp(-r/u) + £y(a, 6, u) exp[-(6 - r)/u]} dw 

(36) 
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in the finite case. 



We do retrieve ip(a, 0) = ip(a, b, 0) = 1 in accordance with (29), due to relations 
(15) and (26). At r = b, Eq. (36) reproduces Eq. (32) owing to (27). 

Expression (35) of the function if) was first given by Minin [12], apart for an 
obvious simplification. It appears also in Heaslet et Warming [14]. Its coun- 
terpart (36) in a finite space seems to be new. 

The end of the calculation of the Q-function is straightforward: in a half-space, 
one obtains from Eqs. (30) and (35) 



where 



Q(a, r) = -—[I - VT^F(a, r)], (37) 

X CL 

' --[1 + F(a,0)-F(a,r)}, (38) 



F(a ' T) - g(a,iA(«)) exp[ - Ma>rl + 2 L wii) eM - T/v)dv - (39) 



The surface value of F follows from Eq. (15) 



F(a,0) = -^=-l, (40) 
V 1 — a 

which justifies Eq. (38). Further, we have 

Q(a,0) = -^=. (41) 

This important result was derived for the first time by Sobolev [10]. It was 
completed by Minin [12], who deduced the first internal expression (37) of 
the Q-function from his calculation of the resolvent function in a semi-infinite 
atmosphere. 

In a finite medium, we introduce the following functions F + and F_: 

F±(a,b,r) = 

R(a){l - £ x [a, b, l/k(a)\ ± £ Y [a, b, l/k(a)]}{exp[-k(a)r] ± exp[-A;(a)(fo - r)]} 

+ o / 9(a,v)[l - £x(a,b,v) ± ^ y (a, b, v)]{exp(-r/v) ±exp[-(6 - r)/v]}dv, 
2 Jo 

(42) 

which yield F again as b — > +oo. Adding and subtracting Eqs. (26) and (27), 
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one obtains 

1 - |(a ± A))M) = (1 - <0[1 + F T (a, b, 0)], (43) 
so that, from Eq. (23) 

(l-o)[l + F + (a,6,0)][l + F_(a,6,0)] = 1. (44) 

The surface values of the functions F± can be deduced from Eqs. (43) and 
(24)-(25): 

F + (a, b, 0) = F + (a, b, b) = ) - 1, (45) 

(1 — a)X (a, b, +oo) 

-F_(a, 6, 0) = F_(a, 6, 6) = 1 - X(a, 6, +oo). (46) 

The counterpart of Eqs. (37)-(38) in a finite slab is, from Eqs. (31), (32) and 
(36) 

Q(a, b, r) = -J— [1 - (1 - a)X(a, b, +oo)F+(a, b, r)], (47) 
= X(a, b, +oo)[l + F + (a, b, 0) - F+(a, b, r)], (48) 

the second equation resulting from Eq. (45). We have furthermore 

Q(a, b, 0) = Q(a, b, b) = X(a, b, +oo). (49) 

From Eqs. (24)-(25) and (43), X(a,b, +oo) is given by 

X(a, 6 , + oo)=l + F_( a ,M)= (1 _ a)[1+ V +(ttM)1 . (50) 

Since X(a, +oo, +oo) = 1/y/l - a, Eqs. (47)-(50) yield Eqs. (37)-(41) again 
when b — > +oo. 

The surface expression (49) - with X(a,b, +oo) as given by Eqs. (24)-(25) 
- was derived for the first time by Sobolev [11]. The first internal expression 
(47) is new, and the second (48) was given by Danielian [15], using the method 
summarized in the next section. 



4 Sectionnd method for solving the problem (5) 



In view of solving the slab albedo problem, Ambartsumian [16] introduced the 
auxiliary function B — B(a, b, r, u) solution to the following integral equation: 
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B(a,b,r,u) = exp(-r/«) + ^ / E 1 {\r - r'\)B(a, b, r', u) dr'. (51) 

2 Jo 

(l/47r)B(a,b,T,u) is the mean intensity of the total (direct + diffuse) field as 
produced at depth r by parallel rays incident on the top surface r = at an 
angle arccos u to the inward normal with unit flux per unit area normal to the 
rays. It is well known [4] that Eq. (51) has a sense and a unique solution for 
any u > 0, thus for u — > +oo. It coincides then with Eq. (5), so that one can 
write 

Q(a,b,r) = B(a,b,T,+oo). (52) 

A recent review on the analytical calculation of the _B-function is in [9], whose 
Eqs. (64) and (66) lead to 

B(a, t, +oo) = H(a, +oo) lim [ui](a,T,u)] (53) 

in a semi-infinite medium, and to 

B(a, b, r, +oo) = X(a, b, +oo) lim [urjxia, b, r, u)] 

— Y(a, b, +oo) lim {wnyia, b, r, u)] (54) 

u— >+oo 

in a finite slab. The limits of the functions urj, ur] X and urjy follow immediately 
from the expressions (69), (74) and (73) of the functions 77, i]x and t]y in [9]. 
Substituting these limits into Eqs. (53)-(54) and using the relations (10) and 
(24)- (25), we retrieve Eqs. (37) or (47) depending on whether the atmosphere 
is semi-infinite or finite. 

Danielian [15] proceeded that way to calculate the Q-function (he denotes by 
N): his relation (22), whose left-hand side contains a sign error, coincides with 
our expression (48) of the Q-function. 



5 Further developments 

Computing the Q-function with the help of Eqs. (37)- (38) in a half-space, or 
Eqs. (47), (48), (50) in a finite slab, is not easy for at least two reasons: 1) these 
expressions contains the indeterminate forms Oxoooroo — 00 as a — > 1, since 
the functions F and F + diverge when a — > 1, while F_ and X remain bounded, 
2) in the finite case, the calculation of the functions F± from their definition 
(42) is time consuming. It requires the prior calculation of the functions £x 
and £y as defined by Eqs. (20)-(21), which presupposes the knowledge of the 
X- and F-functions over [0,1]. 
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To get round the first difficulty, we adopt the expressions (38) and (48) of the 
Q-function in preference to the expressions (37) and (47), and substitute the 
definitions (39) and (42) of the functions F and F± respectively. We obtain 

Q(a,r) = — ^[l + r/(a,r)] (55) 
V 1 — a 

in a half-space, where 

k(a)R(a) „ a ^(a,u) . . du , . 

/(fl ' T) = %lA(a)] 7 M(a)r] + 2io Wia^v) 1 ^ T ^~ ^ 

and 

Q(a, 6, r) = X(a, b, +oo) [1 + r(6 - r)/(a, 6, r)] (57) 
in a finite slab, where X(a, b, +oo) is given by the first Eq. (50), viz. 

X(a,b, +oo) = 1 

+ b^k(a)R(a) [1 - &(a, 6, l/fc(a)) - ^(a, 6, 1/A;(a))] 7* [1, fc(a)r] 

+ -/ </(a,T;)[l-Ma,6,T;)-Ma,6,«)]7*(l,6/«)— \, (58) 
2, Jo f J 

and 

f(a, b, r) = A; 2 (a)i?(a) {1 - &[a, b, l/k(a)\ + fr[a, b, l/k(a)]} 

xf[U(a)r]f[l,%)(6-r)] 

dv 



+ °- f 1 g(a, v) [1 - (a, 6, u) + £y («, 6, v)] 7*(1, r/u) 7*[1, (6 - r)/v\- 
I Jo 

(59) 



t> 2 



We have introduced the function 

7*(M) = " (l " e~ x ) ' (6°) 

which belongs to a family of classical special functions derived from the in- 
complete gamma function: see e.g. Abramowitz & Stegun [17], Section 6.5. We 
note that the values of the function 7* are in the range [0, 1] for x > 0, and 
that 7*(1,0) = 1 and 7*(l,x) ~l/i^0asa;^ +00. A useful expansion of 
7*(l,x) for small values of x is 



x n 



(61) 



The expression (55)- (56) of the Q-function in a semi-infinite atmosphere is our 
final one, since the if -function is easily computable over [0, 1] (Section 6). On 
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the other hand, the formulae (57)-(59) involve the functions £x and £y, which 
are difficult to compute over [0, 1]. The problem of calculating these functions 
is the subject of a wide literature since the fifties [7,18,19,20,21,22]. Recently, 
Rutily et al. [23] have introduced two auxiliary functions ( + and depending 
on parameters a, b and variable u G K, which allow to calculate the functions 
X, Y, £x and £y in one step. Their calculation algorithm is summarized in 
the next section. Formulas expressing the functions £x and £y are, for u > 

1 - &(a, b, u) = -J_I(C + + C_)(a, 6, «), (62) 

M«, 6, «) = 777^(C+ - C-)(a, 6, «). (63) 
H (a, m) 2 



Consequently 



l-£,x(a,b,u)± £ Y (a, 6, u) = — -(± (a, 6, u) , (64) 

tl (a, w) 



and the expressions (58), (59) of X(a, b, +oo) and f(a, b, r) read now 
X(a, b, +oo) = 1 + b {^^]C-[a, 6, !/*(«)] 7*[1, fc(a)6] 



2 /o 1 ^-(.M) 7 UV<}, (65) 



/(fl ' 6 ' r) = g[a,l/fc(a)] fc(fl)C+[a ' 6 ' 1/fc(fl)] 

xf[l,fc(a)r]f[U(a)(6-r)] 

+ 1 1 l^ C+(a ' 6 ' " } 7 * (1 ' r/u) 7 * [1 ' {b ~ T)M ^- (66) 



When b — > +oo, these formulae yield X(a, +oo, +oo) = l/y/l — a and 

/(a, +oo,r) = /(a, r) [given by (56)] as it should, since C±( a ;^) M ) — *■ 1 from 

Eqs. (70)-(74) below, and 67* (1, aft) -> l/a from Eq. (60) (a > 0). 



6 Numerical calculations 



The numerical evaluation of the Q-function with the help of relations (55)- 
(56) and (57), (65), (66) requires the previous calculation of four auxiliary 
functions: k(a), H(a,u), ( + (a,b,u) and £_(a, 6, u) for u G [0,1] and at u — 
l/k(a). The functions R(a), g(a,v) and 7*(l,x) are defined by Eqs. (13), (14) 
and (60). 
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6.1 Coefficient k(a) 



We used the analytical expression of k(a) first given by Carlstedt and Mul- 
likin [19], viz. 



k(a) = VT 



a exp 



9(a, v) dv/v 



where 



9 (a, v ) = — arctan 

7T 



it av 



2T(a,v) 



(0<v< 1). 



(67) 
(68) 



Continuous values on [0, ir] of the arctan function are used in the above defi- 
nition of the ^-function, i.e., the branch is not the principal one. The resulting 
function arctan(x/y) is usually denoted by ATAN2(x, y). With this choice, 
the function v — > 9(a,v) is continuous from [0, 1[ to [0, 1[, although T(a,v) 
vanishes once in the interval [0, 1[. 



6.2 Function H (a, u) 



There are many analytical expressions of this function in the literature, in- 
cluding the following one [18] 



H(a, u) 



we have selected. 



l + u 



1 + u k(a 



exp 



u 



l l 9{a,v)- 
Jo v 



dv 



(v + u) 



(u > 0) (69) 



6. 3 Functions (± (a, b, u) 



The definition and the method of calculation of these functions are explained 
in two papers in preparation [24,25]. Here, we just give the useful formulas, 
with no justification. 



For any u > 0, one has 

£±(a, b, u) = p±(a, b, u) ± M±(a, b)a±(a, b, u), 

where p± and a± are solution to the integral equations 

i i \ ■, i a f 1 9( a i v ) / , / x / , \ dv 
p±{a,b,u) = l±-u ^ ^ exp(-b/v)p±(a,b,v) 



V + u 



(70) 



(71) 
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2k(a)u a f 1 g(a,v) . . . , . dt> . . 

<r± a, 6, u) = ± -u / r exp(-b/v)a± a, 6, v) — — , 72 

1 + k(a)u 2 7o H 2 (a,v) v + u 

which are of Fredholm type on [0, 1]. Once they have been solved on this inter- 
val, these equations define the functions u — > p±(a,b,u) and u — > a±(a,b,u) 
on [0, +oo[, in particular at u — l/k(a). This yields the coefficient M±(a,b) 
appearing in Eq. (70), since 

M ± (a, b) = ^^^fS „ (73) 
±K ' It q(a,b)a±[a,b,l/k(a)Y K } 



where 



<74) 



Calculating the functions (± thus requires the numerical solution of Eqs. (71) 
and (72). Actually, it is proved in Ref. [24] that the solution of Eqs. (72) can be 
reduced to that of Eqs. (71), i.e., the functions a± are analytically expressible 
in terms of the functions p±. It follows that Eqs. (71) are the only ones to be 
solved numerically. Details are not given here, since there is no inconvenience, 
at the level of accuracy, to solve both sets of equations numerically. Indeed, the 
free term and kernel of Eqs. (71)- (72) are regular over [0, 1] and [0, 1] x [0, 1] 
respectively. Their solution p± and a± are regular and smooth over [0, 1], i.e., 
not only continuous, but derivable everywhere in [0, 1], including at 0. 

We note that our algorithm for calculating the functions 1 — £ x ± £y with 
the help of Eqs. (64) and (70)-(74) was already introduced, apart from a few 
details, by Rutily [20] and Danielian [21,22]. 

Every integral in the above expressions has been calculated using the routines 
of the NAG (Numerical Algorithms Group) Fortran Library (D01AHF and 
D01AJF). We also took the routine for solving the Fredholm integral equa- 
tions (71)-(72) from this library (D05ABF). We have collected in [25] some 
arguments showing that these integral equations are currently solved with an 
accuracy better than 10~ n for any value of a and b. 

For given functions H(a,v) and (±(a,b,v), the numerical evaluation of the 
integrals in Eqs. (56), (65) and (66) is difficult in two circumstances: (?) in 
a weakly scattering medium (a close to 0), due to the rapid variation of the 
function v — > g(a, v) in a neighbourhood of v — 1, (ii) when r tends to on 
the right or to b on the left, because of the sharp maximum of the integrand 
close tot> = Tort> = 6 — r (resulting from the presence of the exponentials). 
These difficulties can be overcome thanks to two changes of variable. The first 
one (i) by choosing t — [1 — ln(l — v )] _1 as a new variable, the regularizing 
effect of this transformation being explained in details in the Section 8.1.3 
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of [20]. The second difficulty (ii) is tractable by making the change of variable 
x = [1 — (lnf) -1 ] -1 , since it substantially reduces the strong variations of 
the function v — > exp(— t/v) in a neighbourhood of v — when r — > + 
[26]. When r — > &~, we come down to the preceding case by remarking that 
f(a,b,r) = f(a,b,b-r). 

With these precautions, the integrals appearing in Eqs. (56), (65) and (66) 
can be calculated with an accuracy better than 1CT 11 using the routines of the 
NAG library. Finally, we may safely conjecture that the problem (5) has been 
solved with an accuracy better than 10~ 10 , this estimation being corroborated 
by the curve of residuals of Fig. 1. Residuals are the relative difference between 
the left- and right-hand sides of Eq. (5), calculated for any value of r in [0, b}. 



10" 




10" 5 10" 4 10" 3 10" 2 10" 1 10° 10 1 10 2 10 3 

X 



Fig. 1. Residuals of Eq. (5) as a function of r, showing the accuracy of our solution 
for 1 — a = 10~ 4 and b = 2000. These residuals are symmetric about the r-mid-plane. 



6.4 The case a — > 1 



Until now, the conservative case a = 1 has been excluded, since most of the 
preceding results lead to indeterminate forms when a — 1. This situation 
generates important roundoff errors when a — > 1, which are incompatible with 
the required accuracy. To remove this difficulty, it is appropriate to introduce 
asymptotic expansions of the functions k(a), H(a,u), ... for a — > 1. 

The following expansion of k(a) is given by Case et al. [8]: 



k(a) = y/3(l - a) 



I (1 _ a ) _ HliL!(i _ a f - hi - a f 
5 V ; 5 2 x 7 V ; 5 sV ; 



+ 



2 x 83 



5 3 x 7 2 x 11 



:i-a) 4 + 



(a-1). (75) 
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The calculation of the if-function from its expression (69) does not raise any 
difficulty when a — > 1, except at u = l/k(a). Actually, Eq. (15) clarifies the 
behaviour of the ratio k(a)R(a)/ H[a, l/k(a)} appearing in Eqs. (56), (65) and 
(66), as a — > 1. We have in particular 

k(a)R(a) r- _ 

V3 as a -> 1. (76) 



H[a,l/k(a)] 



As a result, the /-function defined by Eq. (56) tends to /(l,r) = \/3 + 
>/3?(t) — 1 /r as a — > 1, where 



r r 1 g(l,v) dv 

I JO ti (1, f ) f 



(77) 



is Hopf's function [1, p. 139]. This function increases from q(0) = l/y/3 to 
q(oo) = 0.71044608959876; it is tabulated e.g. in [1, p. 141]. 

We conclude that in a semi-infinite atmosphere 

Q(a,r)~ ^3/(1 -a) [r + q(r)} (a - 1). (78) 

In the finite case, there is no difficulty to compute the functions p± and a± from 
Eqs. (71)-(72) as a — > 1. On the other hand, the coefficient M + (o, 6) diverges 
like l/k(a) as a — > 1, while M_(a, 6) remains bounded. Since cr + (a, 6, -u) oc k(a) 
for -u l/k(a), (+{a, b, u) as defined by (70) is bounded when u ^ l/k(a), and 
diverges like l/k(a) for u = l/k(a). Indeed, putting u = l/k(a) and a — > 1 
into Eq. (64) yields 

/3 

fc(a)C+[a,6,l/fc(a)]~^-A)(l,6) as a -> 1, (79) 

where /3 (1, &) is the value at a = 1 of the /^-coefficient defined by the second 
Eq. (22). On the other hand, (_(a,b,u) remains bounded as a — > 1 whatever 
u > 0. 

It follows from Eqs. (23)-(25) that 

x(1 - 6 - +oo) =A(b)- <80) 

Moreover, taking Eqs. (76) and (79) into account, Eq. (66) reads at a = 1 
/(l,6,r) = ^A)(l,6) 



+ 2 



5 i 1 lfe^ C+(1 ' 6 ' V) 7 * (1 ' ^ 7 * [1 ' {h ~ T)/V] ^- (81) 
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We conclude that the Q- function is bounded as a — > 1, with expression given 
by Eqs. (57), (80) and (81) at a = 1, viz. 

«lAr) = 5r(*-r) + s ^{l + r(»-r) 



x i jf ' f^)C + (l, b, ») T*(l, 7» 7 *[1, (* - 'OH^J ■ (82) 



This is an alternative form of Danielian's expression for a = 1 [15, Eq. (23)]. 
The coefficient /3d(1, 6) can be accurately calculated using the following for- 
mula, given here with no proof: 

M1, b) = 73 [6/2 + g(cx))]p_,o(l,6)+P-,i(l,6)' (83) 
where we have introduced the coefficients 

P±>, 6) = 5n,o ± % ^ -^rK exp(-b/v)p±(a, b, v)v n dv (n > 0) (84) 
2 Jo H z {a, v) 

(S n fi = Kronecker symbol). 



7 Tables of the function S — (1 - a)Q 



Numerical data concerning the Q-function are few. We are aware of some 
values published by Sobolev and Minin [27], and of far more detailed tables 
given by Buell et al. [28] using the invariant embedding method. In view of 
checking their five-digit tables, we have computed the Q-function for those 
values of a and b they have selected (Table 1). We got a poor agreement for 
small values of r, most of time 3 digits in agreement only. 

It is clear that the function Q is appropriate for solving analytically the prob- 
lem (1) in a homogeneous and isothermal atmosphere, the solution being given 
by (4). From a numerical point of view, the one drawback of this function is 
that it diverges like — a when a — > 1, b — > +oo and r — > 0. That's why 

we chose to tabulate the function S — (1 — a)Q solution to Eq. (6). It is the 
function that is generally calculated as a benchmark solution to the problem 

(!)• 

The surface values of the S'-function follow from Eqs. (41) and (49), which 
lead to the famous a/1 — a- law in a half-space 

5(0,0) = Vl - a, (85) 
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Table 1 

Values of Q(a,b,r) for a, b and r taken from Tables 1-4 of Buell et al. [28], where 
they are denoted as A, x and t respectively. 



a, b 

7 


T 


Q(a, b, t) 


a, b 

7 


r 


Q(a, b, t) 


0.2, 1 





1.0964212391 


0.9, 4 





2.8413250008 




0.005 


1.0994952541 




0.5 


4.5154025141 




0.U1 


1.1U1o1o0o69 




1 


5.3/ loOUboOo 




n k 
U.o 


i.io4zooioyy 




o 
L 


n OQ Q/i /lono 

0.y<504o44yUZ 


0.2, 5 





1.1178464365 


0.9, 10 





3.1478960008 




0.5 


1.2013059145 




0.5 


5.1418113055 




1 


1.2259902468 




3 


8.6081352918 




1.5 


1.2368227109 




5 


9.1352617263 




Z 


l.z410oobU to 


1, 1 





2.0673668738 




K 

z.o 






0.1 


2.4025482660 


0.6, 1 





1.3962004895 




0.2 


2.5919961827 




0.1 


1.5178886319 




0.3 


2.7148715436 




0.2 


1.5816080857 




0.4 


2.7852328124 




0.3 


1.6213474973 




0.5 


2.8082221832 




0.4 


1.6435798894 


1, 10 





9.8907816365 




0.5 


1.6507643063 




1 


25.586442568 


0.6, 5 





1.5769290982 




2 


36.264849535 




1.5 


2.3199770925 




3 


43.801265888 




2 


2.3672862917 




4 


48.310055291 




2.5 


2.3816047458 




5 


49.811793274 



and to 

S(a, b, 0) = S(a, b, b) = (1 - a)X(a, 6, +oo) (86) 

in a finite slab. Table 2 contains surface values S(a, b, 0) for a wide range of 
parameters a and b, which provides a good test for codes solving the transfer 
equation in a finite slab. 

Within the atmosphere, the S'-function is calculated by multiplying both sides 
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Table 2 

Surface values S(a, b, 0) for a wide range of parameters a and b. 



b 


1 _ a = 10~ 10 


1 _ a = 10~ 8 


l _ a = 10~ 6 


2 


2.9549172423E-10 


2.9549171567E-08 


2.9549086009E-06 


20 


1.8551036585E-09 


1.8551015456E-07 


1.8548902847E-05 


200 


1.7443543255E-08 


1.7441791865E-06 


1.7268758494E-04 


2000 


1.7331077826E-07 


1.7161299590E-05 


9.3944242810E-04 


2 x 10 4 


1.7150550856E-06 


9.3931230259E-05 


1.0000000000E-03 


2 x 10 6 


1.0000000000E-05 


1.0000000000E-04 


1.0000000000E-03 


2 x 10 8 


1.0000000000E-05 


1.0000000000E-04 


1.0000000000E-03 


b 


1 _ a = 10" 4 


1 - a = 0.01 


a = 0.9 


2 


2.9540533187E-04 


2.8714541813E-02 


2.3161893202E-01 


20 


1.8340522066E-03 


9.5164307616E-02 


3.1622023495E-01 


200 


9.4072188464E-03 


1.0000000000E-01 


3.1622776602E-01 


2000 


1.0000000000E-02 


l.OOOOOOOOOOE-Ol 


3.1622776602E-01 


2 x 10 4 


1.0000000000E-02 


l.OOOOOOOOOOE-Ol 


3.1622776602E-01 


2 x 10 6 


1.0000000000E-02 


1.0000000000E-01 


3.1622776602E-01 


2 x 10 s 


1.0000000000E-02 


1.0000000000E-01 


3.1622776602E-01 


6 


a = 0.5 


a = 0.1 


a = 10~ 3 


2 


6.8787981968E-01 


9.4657359383E-01 


9.9948108678E-01 


20 


7.0710678074E-01 


9.4868329804E-01 


9.9949987494E-01 


200 


7.0710678119E-01 


9.4868329805E-01 


9.9949987494E-01 


2000 


7.0710678119E-01 


9.4868329805E-01 


9.9949987494E-01 


2 x 10 4 


7.0710678119E-01 


9.4868329805E-01 


9.9949987494E-01 


2 x 10 6 


7.0710678119E-01 


9.4868329805E-01 


9.9949987494E-01 


2 x 10 8 


7.0710678119E-01 


9.4868329805E-01 


9.9949987494E-01 



of Eqs. (55) and (57) by l - a. We obtain 

S(a, t) = y/l — a [l + r/(a, r)] 
19 



in a semi-infinite atmosphere, and 

S(a, b, r) = (1 - a)X(a, b, +00) [1 + r(6 - r)f(a, b, r)] (88) 
in a finite slab. 

The S'-function has been calculated for the following couples (a, b): (0.5, 2), 
(0.99, 20), (1 - 10~ 4 , 2 x 10 3 ) and (1 - 10~ 8 , 2 x 10 8 ) (Table 3). The first 
value represents an "easy" case, the second, the third and the fourth ones are 
typical of a stellar continuum, an "average" line and a strong line respectively. 
All data are given with eleven figures, which we hope to be significant in view 
of the residuals of Fig. 1. 



8 Conclusion 



Many codes for solving the radiative transfer equation in a stellar atmosphere 
make use of the idealized problems (5) or (6) as a computational test. The 
simplest test is provided by the surface value of the solution S to (6), which is 
close to VI — a m a slab of great optical thickness. This well-known result does 
not turn a rich literature on the subject to best account, which has motivated 
the writing of the present paper. 

From an analytical point of view, we have retrieved the solution to Eq. (5) 
in a semi-infinite atmosphere, which can be attributed to Sobolev [10] and 
Minin [12]. In a finite atmosphere, we have completed the solution given by 
Rogovtsov and Samson [13] and Danielian [15]. Our approach is classical up 
to and including Section 4, as it was developed essentially in the countries of 
the ex-USSR since the end of the fifties. 

In Section 5, the solution to Eq. (5) has been expressed in terms of two non 
classical auxiliary functions, denoted by ( + and which are supposed to 
replace the functions X, Y, £x and £y. Actually, there are some evidences 
that these four functions are not the best intermediate step for solving the 
radiative transfer equation within a finite slab. In a paper in preparation [29], 
we will describe a more direct solution to Eq. (5), which does not involve the 
auxiliary functions 0, ip, B, etc. 

Finally, accurate tables of the surface and internal values of the function S = 
(1 — a)Q are given in this article. We have already taken advantage of these 
tables in a recent article testing the accuracy of the ALI code, widely used in 
stellar atmospheres modelling [30]. 
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Table 3 

S(a, b, t) as a function of r for several couples (a, b). 



a = 0.5 


b ■- 


= 2 


a = 0.99, b = 20 


r 




S(a,b, t) 


r 


S(a,b, t) 







6.8787981968E-01 





9.5164307616E-02 


10~ 4 




6.8805116359E-01 


10~ 4 


9.5218140055E-02 


5 x 10" 


-4 


6.8859830701E-01 


10~ 3 


9.5594616111E-02 


10~ 3 




6.8919789714E-01 


n ni 

VJ.VJ ± 


Q 8408033401 F.-09 


X 1U 


-3 


D.youy4u4 / oyiij-ui 


0.1 


1.1789213897E-01 


0.01 




6.9713223336E-01 


0.5 


1.8112719045E-01 


0.05 




7.2064553226E-01 


1 


2.4654055042E-01 


0.1 




7.4199582863E-01 


2 


3.5613573437E-01 


0.5 




8.2863857712E-01 


5 


5.7511174272E-01 


1 




8.5687200196E-01 


10 


6.9556824233E-01 


a = 1 - 


- 10 


" 4 , b = 2000 


a = 1 


- IO" 8 , b = 2 x 10 8 


r 




S(a,b, t) 


r 


S(a,b, t) 







1.0000000000E-02 





1.0000000000E-04 


io~ 4 




1.0005867102E-02 


io- 4 


1.0005884722E-04 


10~ 3 




1.0047210816E-02 


io- 3 


1.0047386240E-04 


0.01 




1.0359952897E-02 


0.01 


1.0361724723E-04 


0.1 




1.2588036202E-02 


0.1 


1.2607724794E-04 


0.5 




2.0308299268E-02 


1 


2.9416052346E-04 


1 




2.9070093404E-02 


10 


1.8533929890E-03 


5 




9.4251967156E-02 


100 


1.7292311586E-02 


10 




1.6938939655E-01 


1000 


1.5913835240E-01 


100 




8.2524627065E-01 


10 4 


8.2310056337E-01 


1000 




9.9999994061E-01 


10 8 


1.0000000000E-00 
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